clear all
*cap log close
set more off
set seed 603

* appendix version of first stage plot: all outcomes ---------------------------

* setup outcomes --------------------
use "${data}/out/4-main", clear
gen fine100  = fine/100 
label var fine100 "Fine Amount (100s)"
gen efine100 = efine/100 
label var efine100 "Fine Amount (100s)"
label var epoints "DL Points"
cap drop fine 
cap drop efine 
rename fine100 fine 
rename efine100 efine 
qui compress
saveold "${temp}/temp_data", replace

local label1 = "Statutory Fine Amount (100s)"
local label2 = "Paid Fine Amount (100s)"
local label3 = "Statutory DL Points"
local label4 = "Assessed DL Points"

* loop over outcomes ------------------
local k = 1 
foreach D in fine efine points epoints {
	
	* set outcome ---
	use "${temp}/temp_data", clear
	gen D = `D'

	* regression estimate ---
	reghdfe D Z, absorb(totfe) vce(cluster officerid) 
	local b = "{&beta} = `:di %5.3f _b[Z]' (`:di %5.3f _se[Z]')"
	test Z 
	local f = "F = `:di %6.0f r(F)'"

	* residualize for plot ---
	reghdfe D, absorb(totfe) nocons resid 
	predict rD, resid 
	reghdfe Z, absorb(totfe) nocons resid 
	predict rZ, resid 

	* manual binscatter ---------
	binscatter rD rZ, nq(20) nodraw gen(bin)
	gen xbin = _n if _n<=20 
	gen mZ = .
	gen mD = . 

	forval i = 1/20 {
		qui summ rZ if bin == `i'
		qui replace mZ = r(mean) if xbin == `i'
		qui summ rD if bin == `i'
		qui replace mD = r(mean) if xbin == `i'	
	}

	* nonpar first stage -----------
	_pctile rZ, n(100)
	scalar loZ = r(r2)
	scalar hiZ = r(r98)
	lpoly rD rZ if rZ>loZ & rZ<hiZ, nogr ci gen(x y) se(yse) degree(2) bw(0.15) kernel(tri)
	cap gen ylb = y - 1.96*yse
	cap gen yub = y + 1.96*yse

	* build plot --------------------
	#delimit ;
	twoway 
	scatter mD mZ, yaxis(2) msymbol(Oh) mcolor(dknavy) msize(medlarge) ||
	line y x, yaxis(2) lcolor(forest_green) ||
	rarea yub ylb x, yaxis(2) fcolor(navy%30) lwidth(none) lcolor(navy) ||
	hist rZ if rZ>loZ & rZ<hiZ, frac bin(30) gap(20) lcolor(dknavy) fcolor(none)
	graphregion(color(white)) plotregion(lcolor(black) lwidth(medthin)) 
	legend(ring(0) pos(11) order(1 2 3) rows(3) region(lstyle(border))
	lab(1 "Local Mean") lab(2 "Nonparametric Fit") lab(3 "95% CI"))
	xlab(-0.5(0.25)0.5,nogrid) ylab(0(0.05)0.15,nogrid axis(1)) ylab(-0.75(0.25)0.75,nogrid axis(2))
	xtitle("Officer Stringency") 
	ytitle("Fraction of Citations", axis(1)) ytitle("`label`k''", axis(2))
	text(0.055 0.18 "`b'" "`f'", place(e));
	#delimit cr 
	graph export "${out}/apx_iv/fs_`D'.pdf", replace
	local ++k

}

rm "${temp}/temp_data.dta"
* --------------------------------------------------------------------------


* main text version: harsh fine only ---------------------------------------
use "${data}/out/4-main", clear
gen D = harsh 

* regression results ------------
reghdfe D Z, absorb(totfe) vce(cluster officerid)
local b = "{it:{&beta}} = `:di %5.3f _b[Z]' (`:di %5.3f _se[Z]')"
test Z 
local f = "{it:F} = `:di %6.0f r(F)'"

* residualize for plot ------------
reghdfe Z, absorb(totfe) nocons resid 
predict rZ, resid 
reghdfe D, absorb(totfe) nocons resid 
predict rD, resid 

* manual binscatter ----------------
binscatter rD rZ, nq(20) nodraw gen(bin)
gen xbin = _n if _n <= 20
gen mD   = . 
gen mZ   = .

forval i = 1/20 {
	qui summ rZ if bin == `i'
	replace mZ = r(mean) if xbin == `i'
	qui summ rD if bin == `i'
	replace mD = r(mean) if xbin == `i'
}

* nonpar first stage -----------------
_pctile rZ, n(100)
scalar loZ = r(r2)
scalar hiZ = r(r98)
lpoly rD rZ if rZ>loZ & rZ<hiZ, nogr ci gen(x y) se(yse) degree(2) bw(0.15) kernel(tri)
cap g ylb = y - 2.58*yse
cap g yub = y + 2.58*yse

* build plot ---------------------------
#delimit ;
twoway 
scatter mD mZ, yaxis(2) msymbol(Oh) mcolor(dknavy) msize(medlarge) ||
line y x, yaxis(2) lcolor(dkgreen) ||
rarea yub ylb x, yaxis(2) fcolor(dkgreen%30) lwidth(none) ||
hist rZ if rZ>loZ & rZ<hiZ, 
frac bin(30) gap(25) fcolor(none) lcolor(edkblue)
graphregion(color(white)) plotregion(lcolor(black) lwidth(medthin)) 
legend(ring(0) pos(11) order(1 2 3) rows(3) region(lstyle(border))
lab(1 "Local Mean") lab(2 "Nonparametric Fit") lab(3 "99% CI"))
xlab(-0.5(0.25)0.5, nogrid) ylab(0(0.05)0.15,nogrid axis(1)) ylab(-0.75(0.25)0.75,axi(2))
xtitle("Officer Stringency") 
ytitle("Fraction of Citations", axis(1)) ytitle("Pr(Harsh Fine)", axis(2))
text(0.07 0.17 "`b'" "`f'", place(e) size(medlarge));
#delimit cr 
graph export "${out}/main/validity_fs.pdf", replace 
* --------------------------------------------------------------------------


















